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Muon backgrounds at Super-Kamiokande, KamLAND and CHOOZ are calculated using MUSIC. A 
- - ■ modified version of the Gaisser sea level muon distribution and a well-tested Monte Carlo integration 

jir I method are introduced. Average muon energy, flux and rate are tabulated. Plots of average energy 

^^ , and angular distributions are given. Implications for muon tracker design in future experiments are 

I^ ' discussed. 

^ . I. INTRODUCTION 

<: 

ly-s ' Muons are one of the most significant sources of background for underground experiments. An accurate and efficient 

p\j , numerical method to calculate the muon rate and average energy at an underground lab is indispensable for this type 

of research. This work was originally motivated by a need to resolve the question of the average muon energy for 

CN ■ Daya Bay and KamLAND. Since Super-Kamiokande (Super-K) is essentially next to KamLAND and already has 

many publications quoting its muon rates, it easily becomes an ideal source of cross checks. At the same time, the 

need to understand the cosmic background at the far site of Double Chooz also arises. Muon data from the first 

CHOOZ experiment are subsequently made available so that comparison with simulation becomes possible. Given 

^i the diversity of the experimental sites being discussed, some effort is made to present the analysis in a more general 

f^ context. It is hoped that the method presented here will be useful to a larger community. 

^\^ The muon rate can be measured in an experiment by a number of methods. Measurement of muon energy on 

^p I the other hand is quite difficult. Since a measurement made in one site under a certain hill profile is unlikely to be 

(— I . transferable to another site, an economical calculational method is the only practical solution. For these reasons, 

^ i' whenever the average muon energy is needed for the calculation of cosmogenic background rate, accurate Monte 
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JL , Carlo simulation is often the most reasonable alternative. Traditionally there have been some discrepancies in various 
fQ ' reports regarding muon rate and average energy for both Super-K and KamLAND. For example, different values of 
(— I [ muon rates have been reported by different collaborators of Super-K such as 1.88 Hz [l|, 2.2 Hz Q, 2.5 Hz j^j and 
?• 3 Hz ^ . Some of these discrepancies are due to the differences of detector regions or different selection criteria used 

.'^ , in making various cuts. For example, a cut at 1.6 GeV is made to eliminate the muon background in the study 
S^ ■ of the upward throughgoing muons in Reference ;2| . This cut has the effect of lowering the cosmic muon rate. In 
\—t ' Reference |3J, the 2.5 Hz cosmic muon rate quoted is an estimate used to make the spallation cut. Differences in 
cosmic muon rates due to the differences in detector regions will be analyzed by simulation studies later. As far as 
KamLAND is concerned, accurate simulations of the average muon energy, flux and rate are presently needed to aid 
the data analysis process. In addition, the design of muon tracker systems for future experiments depend on detailed 
simulations that can handle complicated topography. In an effort to build a reliable tool for all these needs, this paper 
introduces a complete numerical method from the ground up — beginning with an improved Gaisser sea level muon 
parameterization, showing in detail the logic of the numerical method, making mention of useful numerical tools and 
ending with numerous cross checks with experimental data including those of ground level muons. 

II. A BOTTOM-UP DESIGN 

A. Preliminaries 

The goal of this section is to lay the theoretical foundation for how to incorporate MUSIC with a user-supplied 
sampling algorithm. Details of the implementation of the numerical method outlined in this section can be found in 
Appendix IXI MUSIC is a FORTRAN subroutine that simulates the 3-dimensional transport of muons through a slant 
depth A of a material taking into account energy loss due to ionization, pair production, Brcmsstrahlung and inelastic 
scattering 0,lal- MUSIC is composed of two main parts — (1) the Monte Carlo simulation of muon energy loss and (2) 



the Monte Carlo simulation of angular deviation and lateral displacement. In order to distinguish quantities related 
to initial muons on the surface and the final muons that survive at a certain depth underground, the subscripts /iO 
and fi will be used to denote the two types of muons respectively. 

The testing of the present numerical method involves the comparison of the simulated results against published 
experimental and simulated data. The most convenient item of comparison is the vertical muon intensity IJ^ih) 
versus vertical depth h underneath a flat surface in standard rock because of the abundance of experimental data. In 
order to set the stage for the following discussions, several conventional quantities are defined in the beginning. For 
instance, 6 is defined to be the zenith angle of the line of sight of the muons and (/> is the azimuthal angle of the same, 
measured from the easterly direction in the counter-clockwise sense. Directional muon intensity /^(/i, cos 6) has units 
of cm~^ sr"-'^ s~^. Vertical muon intensity is taken as 

i;:{h) = i^{h, cose ^1), (1) 

which also has the unit of cm^^ sr^^ s^^. Integrated muon intensity is defined as Q 

Mh)= [ I f^iX, cose) dfl, (2) 

where il defines the solid angle coverage and the slant depth X is the distance traveled by the muon in rock. The 
argument X in I^{X, cos0) is generally a function of cos6'. For example, X — h/ cose in the case of a flat surface. 
The unit of integrated muon intensity is cm~^ s~^. 

B. Modified Gaisser Parameterization 

The accuracy of the simulation depends on MUSIC, the parameterization of the surface muon intensity and the 
user-supplied sampling algorithm. A standard atmospheric muon parameterization is given by Gaisser [8| as 
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Muon energy E^^q at the surface is measured in GeV and e is the angle subtended between the incoming cosmic ray 
particle and the normal to the upper atmospheric layer. If the earth were flat, 6* is also the zenith angle on the ground 
surface. Since the earth is not flat, a correlation needs to be made between the zenith angle on the ground surface 
and the angle measured on the upper atmosphere. In order to clarify the distinction between the two definitions of 
angle, a new symbol e* is invented to denote the angle measured on the upper atmosphere as a function of e which is 
assigned specifically to the zenith angle on the ground surface from now on. The calculation of e* will be explained 
later. The symbols A, 7 and Tj. refer to the overall scale factor, power index and the ratio of the prompt muons 
to pions respectively. In the low energy regime, E^q needs to be modified slightly by an energy loss through the 
atmosphere. The symbol Efj_o denotes the energy of muon on top of the atmosphere. The differentials of time dt and 
area dA are omitted from the denominator on the left-hand-side of Eq. ||2Il for the sake of simplicity. The original 
Gaisser parameterization has ^ = 1, 7 = 2.70, Ef^o = E^n and re — 0. For large depth greater than 1-2 km w.e. 
(kilometer water equivalence), the LVD parameterization Q is recommended. In that case, A = 1.84, 7 = 2.77. Since 
this work primarily concerns simulations for relatively shallow depths as in Super-K, KamLAND and CHOOZ, the 
Gaisser parameterization is adequate for the high energy part [E^q > (100/ cos 6*) GeV) of the spectrum. Since there 
are enough low energy muons that survive at shallow depths, rare high energy muons (S^o > 10^ GeV) are omitted 
from the calculations. The valid energy range for the Gaisser parameterization is (100/ cos 6*) < Ef^o < 10^ GcV 
and small angle 6* < 70° where the effect due to the curvature of the earth is negligible. In the low energy limit 
{Efj.0 < (100/ cos 6*) GeV), the Gaisser parameterization is significantly higher than the observed values. The expected 
angular dependence of cos" e with n ~ 2 in this regime must also be taken into account. To satisfy all these 
additional requirements in the small E^^q and large e regimes, the following modifications to Eq. (PJ are suggested for 
(l/cosr) < Ef,o < (100/ cos r) GeV: 

Tc = 10-^ (4) 

A = 2.06 X 10-3 f-^ - 9oV (5) 

Vcos6'* y ' ^ ' 

E^o = E^,o + \ (6) 
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A = l.ll — . (7) 
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It is important to note the term involving cos 9 inside the square root of Eq. Q does not have a star. The LVD 
publications set the upper limit on the ratio of prompt muons to pions to be r^ < 2x 10^'^ at 95% confidence level |9lll0|. 
However x^ of the fits is lower for smaller values of Tc such that the choice of Eq. Q is justified by statistical reason. 
The symbol A in Eq. Q has the interpretation of mean energy loss of muons in the atmosphere. The value 2.06 x 10"'^ 
refers to the stopping power of matter against muons in units of GeV per g/cm^ at -E^o — 50 GeV where the radiative 
effects reach 1% [ll|. The multiplication of this value with the mean muon slant depth in the atmosphere will 
give the mean energy loss of muons in the atmosphere. A commonly quoted value of the atmospheric height hp is 
1000 g/cm^. Reference 12] quotes a specific value of hp — 1030 g/cm^ along with a value of interaction length (the 
average distance between the point where a primary proton enters the atmosphere and the point where a muon is 
produced) Atv = 120 g/cm^. The atmospheric height hp is a function of scale height /iq which in turn is a function of 
temperature. In addition the stopping power used in Eq. jSJl is simply a rough estimate. For these reasons, hp should 
be adjusted to produce the best fit of experimental data in the low energy regime. For the purpose of constructing 
A, the choice of hp = 950 g/cm^ is made. Beside the aforementioned value of interaction length Xn — 120 g/cm^ 
quoted in Reference [l3|, many other values have also been quoted in the literatures, e.g. Xn = 77.6 g/cm^ [I^ and 
Xn = 80 g/cm^ [lj|. Again for the purpose of constructing the fits, a median value of Xn = 90 g/cm^ is chosen for 
A. Putting all these values together, the mean muon slant depth is (950/ cos 6*"^ — 90) g/cm^ such that the final form 
of Eq. jSJ is obtained. For low energy muons, there is a slight difference. A, between the muon energy at ground 
level Ef^Q and the muon energy on top of the atmosphere E^q ■ Since the critical energy (the threshold by which the 
mechanism changes from radiation losses to ionization losses) used in Eq. Q refers to E^q, an adjustment is needed 
for the low energy regime as given in Eq. © . The meaning of Eq. Q is essentially the multiplication of an effective 
factor of 1.1 due to the nuclear enhancement of multiplicity 12] and the probability of muon decay S^. The form 
of Sf^ used in Eq. ^ is similar to that in Reference 12] and can easily be derived as follows: the decay probability 
is related to decay length L = —vjTlnR such that Sf^ — exp{—Xp/L) where Xp = hp/ cos 9* is the slant height of 
the atmosphere, v is muon velocity, 7 = 1/Vl — v'^ is the Lorentz factor, r is the muon lifetime and < i? < 1 is a 
uniformly distributed number (iq . In the muon energy regime E^^q > {1/ cos9*) GeV, the approximation w ~ 1 is 
acceptable. Furthermore, if the choice of i? = X^/Xp is made along with the standard substitutions 7 = Ef^o/m^ and 
hpm^/T ~ 1.04 GeV, the original form of S^ in Reference [l^ is recovered. This derivation reveals that the form of 
S^ in Reference |l3 ] does not incorporate any matter and geomagnetic effects which may be important for low energy 
muons. In the present work, nonlinear effects are taken into consideration by assuming a modification to the decay 
length to achieve the best fits of experimental data such that 
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At this point, Eq. JSJl is purely phenomcnological. There is no simple physical explanation for this change other than 
the fact that it fits the low energy muon data. By replacing L with L in Eq. ((SJ and repeating the derivation of S*^ 
above, Eq. iQ is obtained. 

The modifications in Eqs. Q-Q alone cannot fit the data in the lowest energy range. For £"^0 < 1/ cos 61* GeV, 
the basic form of the parameterization is the same as Eqs. © O with the exception that the substitution 

3£;^o + 7sec6i* 
Et,o -^ — (9) 

is made before -E^o is passed to Eqs. (jSJ-lO. The substitution in Eq. is just another phenomcnological tool to 
achieve good fits with experimental data. 

The value of cos 9* is sometimes calculated using a simple geometrical extrapolation assuming that the altitude 
of first interaction is known a priori. The present work takes a different approach by using a more complicated 
extrapolation method described in Reference !l6| that shows how cos 9* can be extracted from an integral equation 
by equating interaction length X{9) = A'(O). In essence the formula below taken from Reference [16] parameterizes 
the numerical solution of the integral equation: 
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where x = cos9, pi = 0.102573, p2 = -0.068287, pa = 0.958633, pi = 0.0407253 and ps = 0.817285. The terms 
involving cos 9 in Eq. ^ must be replaced by cos 9* for consistency. Eq. ([TJ is protected against division by zero 
because cos 61* > 0.103458 for cos 61 > according to Eq. (|10|l . The modified Gaisser parameterization is based on 
the world data set and hence represent an average of the global sea level muon distribution. The geomagnetic field 



affects only the low energy spectrum, typically below 2 GeV for integrated muon intensity T7] and approximately less 
than 20 GeV for vertical muon intensity [l^- The east- west effect is also shown to be negligible at ground level [I3 
by careful simulations. For the purpose of calculating the muon overburden deep underground, geomagnetic effects 
can be ignored because low energy sea level muons will not survive through rock by default. In essence, the present 
parameterization is composed of the union of 3 segments: Ef^Q > (100/ cos 6'*) GeV (the standard Gaisser formula), 
(l/cosr) < E^,o < (100/cos6l*) GeV (Eqs. IJU-llIIl) and E^,o < (l/cose**) GeV (Eq. ©). Figure [l] illustrates the 
quality of the fits between the modified Gaisser formula and experimental data. The goodness of fit tends to degrade 
only at very large angles {9 > 85°). The worst disagreement between experimental data and the parameterization 
in those cases is about 40%. However the worst case scenario of the 40% disagreement occurs only at low energy 
{Efj^Q < 10 GeV) and a relatively small sector at large angles {9 > 85°) so that the integrated spectrum is dominated 
by the very accurate parts of the parameterization at smaller angles. Finally it should be emphasized that the 
modifications to the low energy part of the standard Gaisser formula outlined in Eqs. Q-Q and © will not have 
any significant impact on the simulations of deep underground experiments. Nevertheless an accurate description 
of the low energy part of the sea level muon distribution is important for calculating the muon overburden for sites 
situated at shallow depths such as the near sites of the Double Chooz and Daya Bay experiments. 

C. Modeling Physical Observables 

Directional intensity at depth h is obtained by integrating the Gaisser parametrization and the muon survival 
probability over the initial muon energy iJ^o at the surface as 

I,{X,cos9\cj^)= r dE,o P{E,o,X, 9\ <j^) '^^^^i^^^^co^ ^^) . (11) 

Jo dE^odO. 

The probability function P(Efj^Q, X, 9* , cj)) defines the survival probability of a muon with initial energy E^^ traversing 
a slant depth X from the zenith angle 9* and the azimuthal angle 0. It is emphasized that the symbol I^ as used in 
this paper has a different meaning than I^ in Reference |3ll | in that the latter refers to a differential muon intensity 
containing a probability distribution function P{E^, -E^tOi ^j ^*i 4') which is related to P{E^q, X, 9*, (f) as per 



P{E,,o, X, 9\ 4>)= I dE^ P{E^, E,,o, X, 9\ <t>). (12) 



In Reference |3l|, 2 x lO"^ muons with energies from 0.1 up to 1000 TeV were propagated through 15 km.w.e.of rock. 
The values of P(Efj^, E^q, X) were stored and then integrated numerically using an equation similar to Eq. (|ll|l to 
obtain energy and angular distributions at any particular depth. This work takes a different approach by evaluating 
P{Efj_o, X, 9* , (f)) in situ in the Monte Carlo integration. This approach requires a smaller number of simulated events 
(typically 5 x 10^) and is more versatile when applied to arbitrary hill profiles when P [E i^^q , X , 9* , (j)) must be re- 
evaluated every time the (x, y)-coordinates are changed. In principle the transport of muons from the surface to 
a point underground and vice versa are equivalent as far as the calculation of energy loss is concerned. The most 
important requirement of the present method is the uniform generation of E^q, 9 and (j) as shown in Appendix^ An 
arbitrary energy dependent observable 0^{Efj) can be estimated as 

(0^(X,cosr,0))= , ,y \, ,, r dE.o '^^'"^^'"'^^'^*^ r dE,0,{E,)P{E,,E,o,X^9\<l^). (13) 
7^ (A, cos f ,(p) Jo aii^o"" Jo 

The bracketed quantity on the left-hand-side of Eq. ((T^ represents the average of O^. The bracket will be dropped 
from now on for the sake of simplicity unless ambiguities arise due to the choice of symbols. With Eqs. (|ll|l and H13|l . 
vertical intensity and average energy are 

^'<"' ^ 7^ r ''"' '""''X^n ^ " r '^' ^' ^<^- ^- "• ^- - "■ «^ ("' 

There are many different ways to implement Eqs. H14I) and (|15|l numerically. Appendix ^ describes an efficient and 
accurate Monte Carlo method. Simulated values of IZiJi) beneath a flat surface are compared against experimental 
and simulated data in Figures [3 The results obtained by using the modified Gaisser parameterization incorporating 



Eqs. 01-0 at low energy agrees with experimental data more closely than those using the standard Gaisser param- 
eterization only. Figure ini shows the consistency between simulated and experimental data at shallow depths. (The 
interpretation of Figure O will be discussed more fully in Appendix 0) The integrated muon intensity and average 
energy are 



dN^o{E,,o,cos9*) 
dE^o dQ 
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Jf^ and (£'p) in Eqs. 1)16(1 and H17|l are functions of the location of the point of sampling underneath a topographic 
profile. The arguments of these functions, namely the coordinates of the point of sampling, are understood and 
therefore not displayed explicitly. The brackets around E^ on the left-hand-side of Eq. (|17|l are dropped in the 
following text whenever the reference to average muon energy is clear from the context of the discussion. Integrated 
intensity can be computed in a similar way as vertical intensity. The only difference is that the depth X is now a 
function of 9 and (j>. In the case of a flat surface, the relation takes on a simple form X = h/ cos 9. In the case of an 
arbitrary hill profile, there is no longer any simple relationship among X, 6 and (f>. 

The key of the present numerical method is uniform generation of integration variables which can be achieved in 
reliable ways through various uniform generation algorithms such as the CERNLIB routine RANLUX. The logic of the 
method is relatively simple so that there is little ambiguity of its correctness. All these observations lead to the 
conclusion that MUSIC is sufficiently accurate over the relevant range of muon energy. Simulated integrated muon 
intensity and average energy are compared against published simulations in Tabled 

III. PREPARING THE CALCULATION 

A. Digital Maps 

The starting point of a muon simulation over an arbitrary hill profile is the digital map of the surrounding topology. 
The accuracy of a digital map directly affects the accuracy of the calculation so that a detailed knowledge of the 
hill profile is important. According to a contour map published in Reference J23i KaniLAND is separated from 
Super-K by 187m and its bearing is N66.6E with respect of Super-K. The top of the Super-K tank and the bottom 
of the KamLAND tank are situated at 350m above sea level. Both sites are almost directly underneath the peak of 
Ikenoyama at 1.35 km. Due to their proximity, both sites have very similar muon energy and flux. However the sizes 
of the two detectors are vastly different so that their muon rates scale accordingly. The digital map of the Ikenoyama 
mountain profile around Super-K is extracted from a code used by M. Nakahata originally to calculate the muon 
background for Supcr-K. The Super-K digital map sets its origin at the location of the detector and parameterizes 
coordinates in terms of (p, 0, h). This particular format does not allow a simple coordinate transformation of the 
origin from Super-K to KamLAND. As a result, the digital map of KamLAND is constructed independently from 
a contour map for this work. In order to guarantee a sufficient solid angle coverage for the simulation, both digital 
maps cover circular areas of radius 3950 m. The topological map of CHOOZ is generated from a 2D contour map 
using a shareware called 3DField p34| . A visualization of the digital map over the Ardennes Mountains is shown in 
Figure^ 3DField has the option of generating an ASCII data file containing the {x, y, z) coordinates of the latticized 
hill profile. One side of the CHOOZ detector is beneath a steep hill so that a large range of the solid angle coverage 
is parameterized by a relatively small set of lattice points. In order to increase the density of lattice points over the 
steep section of the Ardennes hill profile, another digital map is created over a smaller area around the detector. At 
the end, both digital maps are spliced together to form one single digital map so that the entire solid angle coverage 
is represented more evenly. 

B. Detector Geometry 

The calculation of the average muon rate depends on the details of the detector geometry. For Super-K, the 
parameters that define the geometry of the cylindrical tank are Lq = 41.40 m and i?o = 19.65 m. An inner volume 
is defined to eliminate the simulations of very small muon track lengths inside the detector geometry that do not 
intersect the active region of the detector. The choice of the inner volume is not critical for the calculation of muon 
rate inside the outer tank i?* . For the purpose of this work, the inner volume used in the simulation of i?* is also the 
inner detector volume of Super-K whose dimensions are L = 36.20 m and R = 16.90 m. The inner and outer tanks of 



Super-K have almost the same aspect ratios so that the muon rate inside the inner tank i?* can be obtained simply 
by scaling i?* according to the ratio of the physical areas ^o of the two tanks. The geometry of the cylindrical tank at 
KamLAND is defined by Lq = 19.68 m and Rq = 9.50 m. The inner spherical volume of KamLAND for the purpose 
of this simulation is taken to be the area bounded by the buffer region with R — 8.25 m. For the simulation of the 
muon rate inside the KamLAND detector volume, a sphere of radius Rq — 6.50 m is used. In the case of CHOOZ, 
the cylindrical tank has the dimensions of Lq = 5.5 m and Rq = 2.75 m. The inner detector is filled with Gd- loaded 
liquid scintillator and has the shape of a short cylinder with hemispherical end caps. Muon rate inside the Gd-loaded 
region is not simulated in the present work. 

C. Rock Composition 

Chemical composition of the rock affects a MUSIC simulation in that two out of three cross section files need to be 
calculated with specific rock data a priori. Table ^1 gives the chemical composition of the Ikenoyama and Ardennes 
rock. The average atomic number and weight are (Z) = 10.13 and (A) = 20.42 for the Ikenoyama rock and (Z) = 11.8 
and (A) = 24.1 for the Ardennes rock. Hydrogen composition is 2.2% for Ikenoyama and negligible for Ardennes. The 
rock density and the radiative length are p = 2.70 g/cm'^ and A = 25.966 g/cm^ for Ikenoyama and p — 2.81 g/cm^ 
and A = 23.3 g/cm^ for Ardennes respectively. The present simulation for CHOOZ takes the approximate chemical 
composition and the average rock density |35l| as inputs. The Ardennes Mountains has a complicated rock density 
profile with a layer of dense rock (3.1 g/cm^) |36l | on the northeast sector. 

In principle complex geological profiles can be incorporated into the MUSIC simulation by a stratified approach 
in which a simulation is segmented according to regions of different densities, average atomic numbers (Z), average 
atomic masses {A) and radiation lengths. Although the stratified approach is possible, it may not be easily achieved 
in practice. Aside from the computational challenge of simulating a complex geological profile, information of the 
geological profile obtained by geological surveys may not be generated with sufficient details to support a realistic 
simulation in the first place. Fortunately the stratified approach can be avoided in many cases. If (Z) and (A) are 
constant and only density varies with depth, the mean density should give the same average muon energy and flux 
as those generated from stratified densities. Varying densities may affect the profiles of angular distributions as in 
the CHOOZ case shown in Section Hvl Radiation length affects only the lateral displacement, which is not under 
investigation in this paper. Small changes in (Z) and (A) (up to 10%) should not seriously affect the muon flux as long 
as the mean values of all layers are found accurately. This work does not attempt to simulate the detailed geological 
profile of the Ardennes Mountains. It is shown in Section llVl that the simulated results due to the simplification of the 
Ardennes geological profile are consistent with the previous CHOOZ measurements within errors and that simulated 
results of a uniform Ikenoyama mountain profile agree with experimental data. 

IV. RESULTS AND DISCUSSIONS 

A. Average Muon Rate 

The calculation of muon rate depends on the the effective area of the detector. The basic strategy of calculating 
the effective area A is to multiply the physical area Aq with the ensemble average of the inner products of randomly 
generated unit vectors fi pointing from an inner volume and the unit normal vectors fg pointing away from the outer 
surface. In this case, the ensemble average also constrains the pseudo survival probability of muons (P) that will 
be defined more precisely by Eq. IJA14|I in Appendix 1X1 Figure ITHl visualizes how the inner products are done. An 
intuitive way to think about the effective area A is 

^ = ^E^-^- (18) 

i=l 

In the case of a cylindrical detector, the physical area is Aq = ttRq + 2ttRqLq. Similarly Aq = AttRq for a spherical 
detector and so on. If Eq. (|18|l is used, the average muon rate i?^ is simply 

i?M - -^M ^- (19) 

The average muon fiux J^ is always sampled at the center of the detector volume in this work. Although the 
macroscopic strategy defined by Eqs. I|18|l and H19|l gives reasonable results, a microscopic strategy to compute the 



muon rate is considered more accurate, namely 

where dA is a differential area element along the the detector wall and f is a unit vector along the muon line of sight. 
Both are functions of position along the detector wall, cos 9* and cj). Appendix 1X1 gives a numerical implementation 
of Eq. (EOl- 

Table IIIII summarizes the main results in terms of average muon energy, flux and rate. The muon rate in the 
outer tank of Super-K generated by the present method is somewhere between the experimental values published in 
References Q and [J. The 3 Hz muon rate quoted in Reference Q is most likely a rounded figure. It is not clear 
if the 1.88 Hz muon rate in Reference jj] refers to the inner detector volume only. If so, it would agree with the 
simulated result very closely. The muon rate at Kamiokande is usually quoted as 0.44 Hz J33. Since KamLAND 
is slightly larger than Kamiokande, the muon rate should be scaled according the ratio of the physical areas ^0 of 
the two tanks which becomes approximately 0.5 Hz. The unofficial measured rates on the KamLAND outer detector 
and the balloon are 0.75 and .21 Hz respectively. They differ from the simulated results by about 10% and 17% 
respectively. The official measured muon rate in the spherical buffer region of radius Rq = 8.25 m is 0.34 Hz 38] and 
the simulated result is 0.396 Hz (14% difference). The muon flux of 0.4 m^-^s^^ quoted by CHOOZ 36] is smaller 
than the simulated exact result by about 35%, which is attributable to the single digit of precision of the quoted 
rate and the approximated geology used in our simulation. The errors in the simulated results in Table lllll are a 
combination of the systematic error from map-making and the statistical error from the Monte Carlo simulation. The 
systematic error of the mountain profile coming from the calculation of the scale that relates physical distance on the 
contour map to the relative distance on the digital map is taken to be 0.5%. The systematic errors of E^, J^^ and i?^ 
are calculated by varying the slant depth X by 0.5% before passing it to MUSIC. The statistical variance is calculated 
in the usual way by varying the random seed. 

B. Energy and Angular Distributions 



The energy distribution dJ^/dE^^ in Figure [S] is defined by the formula 

dJ, 
/o dE^ 



J^^l dE,^ (21) 



and has the unit of GeV~^ cm~^ s^^. Angular and double differential distributions can also be defined in similar ways. 
Appendix fXl describes numerical implementations for various types of distributions. Figure [S] plots the cosmic muon 
energy distributions at Super-K, KamLAND and CHOOZ. The purpose of the figure is to show the global properties 
of the energy distributions of various experiments. Although the distributions look smooth on the log-log scale, the 
fluctuations in the low energy regime {Ef^ < 1 GeV) will become more apparent on the semi-log scale. Fortunately 
the fluctuation in the low energy part of the spectrum on the log scale is suppressed by the smallness of the Jacobian 
that contains a factor of E^ so that the accuracy of the calculations of the average muon energy E'^ and flux J^ are 
not affected. If the energy distribution of stopping muons is needed, generation of iJ^o and cos 6 according to the 
surface muon distribution is recommended. 

Figures IHHHl illustrate the angular distributions of muons. Experience shows that 5 x 10^ simulated events are 
generally adequate to generate reasonably good quality distributions in most cases. The polar angle in the relevant 
plots is defined to be the zenith angle consistent with Eqs. ^ and |(2Jl. The azimuthal angle (f> is set to zero when the 
final muon travels from east to west. The momentum of the final muon is opposite to the line of sight connecting the 
detector and the muon and is defined by and (p. The only exception to the present definition of 9 and (p is Figure |H| 
because the Super-K digital map uses a different convention. Figure [HI compares the cos6 and (p distributions between 
simulations and experiment at KamLAND. Figure ITUl compares the 9 and (p distributions between simulations and a 
cosmic ray experiment done on the CHOOZ site in 1994. The experiment consists of four 1x1 m^ PRC (Resistive Plate 
Chambers) plates separated from each other by 20 cm. The simulation of the experiment defines a muon event as the 
coincidence of any two of the RFC plates being triggered. The difference between the simulated and the experimental 
9 distributions can be explained by the fact that a significant number of the muons coming from the steep section of 
the Ardennes hill profile cannot be detected by a muon tracker composed of top and bottom horizontal plates only. 
In order to measure the muons coming from large zenith angles, additional trackers are needed on the sides. The 
remaining small differences between the simulated and experimental results and the aberrations presumably arising 
from the variation in geology described in Reference |36| are not simulated in this work. The difference between the 
simulated and experimental dJ/dcp for < cp < 150° in Figure lTUl is consistent with an unpublished result in an internal 
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note of the CHOOZ collaboration. Notwithstanding the lack of detailed treatment of smaller features, a macroscopic 
picture emerges by the way of a qualitative comparison in the performance of two types of muon trackers represented 
by the horizontal plate cosmic ray experiment on the CHOOZ site and the muon veto system of KamLAND. The 
simulated cos 6 and <f> distributions agree well with KamLAND experimental data because the muon tracker system 
of KamLAND has full sensitivity over the entire range of the hemispherical solid angle. The disagreement between 
the exact simulation of the 9 distribution and the experimental data measured by horizontal plates of the CHOOZ 
cosmic ray experiment in Fig. IIUI shows that the contribution of muon flux from the sides cannot be neglected in the 
case of a detector located underneath a hilly topology. The obvious exception to this claim will be the case where a 
detector is situated underneath a flat surface so that the slant depth grows with secO. 

Figures [Tnil2l plot the average muon energy versus 9 and </> for KamLAND and CHOOZ respectively. It is noted 
that the differential flux in Figs. |^ and IIUI tends to vary inversely with the average muon energy E^ per angle in 
Figs. [m and ll2l This anti-correlation is intuitive in that average muon energy generally increases with slant depth 
while muon intensity decreases with slant depth. 

V. CONCLUSION 

The described method integrates MUSIC, a modified Gaisser parameterization of the sea-level muon spectrum, and 
a uniform sampling algorithm for the surface topography. The method is efficient, robust, and portable. Given 
sufficiently accurate geological data, the method is capable in principle of predicting muon rates and mean energies 
within a few percent accuracy for depths less than 2.5 km.w.e., as indicated by the error estimates in Table IhH 
In practice, simulations performed using simplified geology assuming uniform rock composition lie within 10~20% 
of observed rates published by Super-Kamiokande and KamLAND, and within 35% of the published flux at the 
geologically more complex CHOOZ site. Although muon simulations for any arbitrary hill proflle have already been 
done many times by other researchers previously, there are very few complete documents approximating a pedagogical 
introduction to the numerical method itself. Although muon rates can be measured in an experiment, muon energy 
is difficult to measure so that knowledge of the average muon energy depends on simulation. For this reason, the 
reliability of the numerical method is very important. In applications such as the estimation of muon background 
in reactor 9i^ experiments, the method of uniform generation of variables can serve as an additional cross check for 
accuracy. Although the standard Gaisser or LVD parameterizations are generally adequate for the simulations of 
the deep underground experiments, the modified Gaisser parameterization is indispensable for shallow depth muon 
simulations. 
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APPENDIX A: TECHNICAL DETAILS OF THE NUMERICAL METHOD 

The quality of the 3D topological map is crucial for an accurate calculation of muon overburden. It is usually the 
first and the most important step. Care should be taken to remove the disconnected parts of the mountain profile. 
If a ray defined by a specific set of [9, </>) passes through disconnected parts of the mountain geometry resulting in 
several different values of slant depth X, the smallest value of X is used for that particular solid angle. 

The range of the energy in the integral goes in principle from to infinity. It is more advantageous to numerically 
integrate over log_E^o instead of E^q. (In this work, log refers to base 10 logarithm and In to base e.) On the other 
hand, integration over E^o will give essentially the same results. The range of numerical integration over muon energy 
E^Q is labeled by the lower and upper bounds, Ei and £'„ respectively. This work chooses not to change the variable 
in such a way to integrate up to E^q — > cx). More specifically, the natural cut-off point ought to be a sharp drop in the 
muon spectrum which in turn correlates with the knee of the primary proton spectrum between 1000 and 10000 TeV. 



The change in the muon spectrum is 5-10 times lower than that so that a reasonable estimate of the upper limit is 
Eu = 10^ GeV. As a practical consideration, it is more computationally efficient to set Ei not strictly as the rest mass 
of muon rrif^ but the minimum muon energy needed to survive the minimum slant depth of a particular geographic 
profile so that CPU time is not wasted in simulating muons that cannot survive by default. 
After the change of variables from E to log(£'), Eq. Ullfl is transformed as 

/^(X,cosr,</))^lnlO /'"" """ d log E^o P{E,o, X, 9*, cj,)E^„ dN,,o{E^o,co8 6*) ^^^^ 

Jlog El CL-tjuO «5 I 

An integral in the Monte Carlo method |42j can be approximated as 

fix, y)dy, ~ (/(x, y)) ■ (?/2 - yi), (A2) 



lyi 

where {fix, y)) is the average of fix, y) over y. With Eqs. HA1|I and (jA2|) . Eq. ((TT1) can be calculated numerically as 



In 10 (log ^„~log^;) ^ dN^oiE^o^cos 6*) 

N ^ ^°' dE^o,dn 



7^(A,co&6^ ,0) ~ 2^^mo* j^ — ^ • (A^j 



The symbol {i} denotes a subset of simulated events corresponding to surviving muons. Information of X and </> on 
the right-hand-side of Eq. (|A3p are defined as inputs in the simulation and are not shown formally. The probability 
function P(£'po, X, 0* , (f) is not explicitly computed in Eq. (|A3p by design. The simplicity of this algorithm translates 
to saving in memory. Since the probability function is not computed explicitly for each combination of E^^q, 9 and (j), 
it is essential that the generation of these integration variables is uniform so that the probability function is calculated 
implicitly when the sum is divided by N in Eq. (|A3(I . As a test of the accuracy of present method, it will be shown 
later that the uniform generation of integration variables gives exactly the same results as those calculated by the 
Gaussian quadrature method in the case of ground level muons. 
Vertical muon intensity and average energy are easily computed as 

™.,^ ^ lnlO(logK-logi;;0 ^ dN,,oiE^^,,cos9t = l) 

™^.x _ In 10 (log g^- log £^0 .^ ^ dN^j,oiE^,oi,coBe* = l) 

""'^^^ - — i^nfsh) — ^ ''''''''' — dE;:^;jn — ■ (^') 

Vertical muon intensity in standard rock simulated with Eq. IJA4|I is compared against experimental data in Figs. [21 
and 01 In reality, standard rock does not exist and is generally not an exact match of real rock profiles in real 
experiments. When measurements are converted from real rocks to standard rock, there are always some questions 
regarding the accuracy of the conversion schemes. For this reason, Figs. (21 and should not be taken as absolute tests 
of the accuracy of MUSIC and the present integration method but merely a relative point of reference. Despite of the 
question of the accuracy of the standard rock experimental data, it is shown in Fig. that the ratios of calculated 
and experimental vertical intensity scatter symmetrically around unity so that the simulated results are said to agree 
with experiments at large. It should be noted that vertical intensity is merely an approximate test and is not the 
central focus of the present work. Integrated intensity on the other hand is really what is needed for the calculation 
of muon overburden in realistic calculations. Integrated muon intensity and average energy can be implemented in a 
similar way as Eqs. HA4P " HA5|I . 

»lnlO(logg„-log£;;) ^^ rfA^^o(£^o.,cos0*) 

-^^^ N ^^^'' — dE::;;dh — ' ^^^^ 



in 



MOi ' 



mnlOOog^^bg^^ dN^oiE^o,,coser) ,._. 

^- - in^ih) Z. ^M. E,,. ^^^-^-^ , (A7) 

where 51 is the solid angle over the integrated hill profile. The average muon energy E^ can be organized into M = 500 
bins along logE^. The subscript j denotes the j-th bin. The numerical implementation of Eq. (|21|) is 

dJf, QM "^^'^ Ef,o^ dN^oiE^m,coae*) . . ^. 



dE^,j N ^ Ei,j dEf,o,dn 



i=l 
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Information of the survival probability is hidden in N{Ef^j) that gives the number of surviving muons in the j-th bin. 
As a consistency check, 



j=i 



and 



^ ^ lnlO(logg„ - logEi) ^ ^, dJ^(fe) 

""^ - — m;ji^) — h ''^^ ^ ^ 

must agree with those obtained by Eq. (|A6|) and ljA7p . Angular and double differential distributions are constructed 
in similar ways. For instance, the cos 6 distribution can be constructed as 

dJ^ _ n MlnlO(logK„ - log£;i) ^(^^^^ ^ dN^,Q{Ef,oi,coa9*) 



E.O. --^ ^' . (All) 



dcos^j 2 N ~i dEf^oidVL 

The factor 0/2 in Eq. IjAlip gives the proper normalization so that the integration over — 1 < cos 6* < 1 gives the 
correct solid angle 51. Similarly the ((> distribution is given as 

dJ^ _ ^ M\nlQ{\ogEu -log El) ""^^ ^ dN ^qJE^q,, cos 0*) 

d(t)j 2tt N ^^ ''" dEf,o,dn ■ ^ ' 

In the case of Eq. (|A12|I . the normalization factor is il/27r so that the integration around < </> < 27r gives the correct 
solid angle Q. There is a subtlety involving the 9 distribution. Since cos 9 (not 9) is uniformly generated in the present 
method, uniform binning in 9 leads to the wrong distribution. The correct bin width must be inversely proportional 
to cos9 or, more precisely speaking, equals M/{N cosd). The factor of l/cos6' exactly cancels the factor of cos6' of 
the Jacobian so that the 9 distribution is 

dJ^ _ n MlnlO(logg., - logg;) "^^ ^ dN^o{E^o,,cos9*) 

d0, " TT N ^ ^"^ dE^o, dn • ^ ^^ 

Double differential distributions of various kinds can be constructed using the same logic. 

Exact calculation of the slant depth X is generally impossible for any arbitrary latticized hill profile. Fortunately 
simulated energy loss by MUSIC is not very sensitive to small changes in X so that an approximation can be made. 
Figureini illustrates the binning strategy of X in the 9—(j) space. The idea is to partition the solid angle into regions of 
nearest neighbors. Each region has the same value of X. Evenly generated 9 and </> pick out an approximate value of 
X in the corresponding nearest neighborhood. There is a certain amount of computation overhead in pre-processing 
the partitions. When the number of simulated events is sufficiently large, the overhead of partitioning is still more 
cost-effective than a real-time search per event. Due to the irregularity of the hill profile, any given differential solid 
angle in the upper hemisphere may traverse disconnected parts of the hill profile. For this reason, the code must 
incorporate a mechanism to pick out the appropriate slant depth X. It can be easily implemented by simply keeping 
only the minimum value of X for any given grid in a latticized 9 — (j) space. On ground level, muons do not need 
to be propagated by MUSIC so that Eqs. (|16|l and (|17|l can be integrated by Gaussian quadrature and by setting 
P{E fj^o, X,9* ,(f>) — 1 in the integrand. Figures [l5l 1 1 61 and Table Hvl show that results generated by the present Monte 
Carlo method agree exactly with those calculated by the Gaussian quadrature method. 

Assuming that pairs of cos0, </> are uniformly generated N times and that the generation is truly random, the 
hemispherical fl would have been partitioned into N segments uniformly. In that case dA ■f = Aof- tq/N so that the 
integral of the differential projected surface areas is simply the sum of the segments corresponding to the surviving 
muons. In other words, for any given f, J dA. ■ f ~ Aq (P) f ■ tq where (P) is the pseudo survival probability of 
muons. A real survival probability of muons can only be computed by generating muons according to the sea level 
niuon distribution and by propagating them through rock. In the present method, muons are generated uniformly in 
E^Q , cos 9 and (p so that (P) does not have any natural meaning other than the ratio of the surviving muons to the 
generated muons. Naively (P) may be set to the ratio of the number of surviving muons n to the generated muons 
TV according to this definition. However a more careful look reveals that n increases as the muon energy threshold Ei 
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is raised in the simulation. The reason is simply that more highly energetic muons generally have a better chance of 
surviving through rock. For this reason, a proper definition of (P) must be independent of Ei and is founded to be 

W-^ ,'°'f-',"''f' - (A14) 



N log Eu — log m 



M 



(P) defined in Eq. (jA14p is effectively independent of Ei because n varies inversely with log Ei. On the other 
hand, n decreases when E^ decreases because less highly energetic muons are generally less likely to survive through 
rock. Unfortunately there is no simple way to rescale (P) in this case. For an accurate calculation of R^{h), it is 
recommended that E^ is kept at 10^ GeV. With {P) in place, Eq. H2U|I can be implemented as 



Ao (P) rjln 10 {\ogEu - log^O y- ~ ^_^ _ dN^^iE^a,, cos 9*) 

{i} 



n aq {r-)\imw [.log i^u - log J^i j .^ , , aiMf,o(i^f,Q.„cosu- ) 

^' " N ^ '' • ''' ^-"^ — dE;:;;;dh — ■ ^^^^^ 



The dot products fi ■ foi are generated randomly inside the entire detector volume and not just at the center. This 
strategy renders a fairer sampling of the detector geometry. On the other hand, the hill profile is defined with respect 
to the origin which is normally set at the center of the detector because the generation of slant depth X is not easily 
managed when the origin moves. Since MUSIC is relatively insensitive to small change in X and the size of the detector 
is generally small compared to the mountain profile, generation of X from the center of the detector is adequate. 
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FIG. 1: Fits of the modified Gaisser parameterization to experimental data in the low energy regime between 6 — and 
9 = 87°. The experimental data are taken from References [g, |20( I-|3C | . The modified Gaisser parameterization is given by 
Eqs. ©-Q and lP- lfTn)l . 
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FIG. 2: Average vertical muon intensity IJ^{h) versus vertical depth h beneath a flat surface in standard rock. Experimental 
and simulated data are labeled by EX and MC respectively. The experimental data of the flat surface overburden are taken 
from References [M |32| and the simulated data by Kudryavtsev et al. from Reference '|3l|. The number of simulated events 
per data point in this figure is A'^ = 10^ . The set of simulated data labeled as "Present work" is generated from the modified 
Gaisser parameterization of the surface muon intensity in Eqs. 0-0. 
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FIG. 3: Ratio of simulated vertical muon intensity J^(/i)mc over the experimental vertical muon intensity I^{h)EX versus 
shallow depth h beneath a flat surface in standard rock. The experimental data of the flat surface overburden are taken from 
References [3,|32|. The number of simulated events per data point in this figure is A'' = 10^. The simulated data are generated 
from the modified Gaisser parameterization incorporating Eqs. Q-Q. 
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FIG. 4: A visualization of the 3D topological profile of the Ardennes Mountains over the CHOOZ site generated by SDField 
from a 2D contour map. 
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FIG. 5: Integrated muon intensity distribution at Super-K, KamLAND and CHOOZ. The number of energy bins is M — 500. 
The total number of simulated events is 5 x 10^. The average muon energies of the three sites are quoted in the legend. 
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FIG. 6; Angular distribution of final muons at Super-K. The total number of simulated events is 5 x fO®. The Super-K digital 
map defines = to be along the northerly axis and the sense of rotation to be clockwise. 
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FIG. 7: Angular distribution of final muons at KamLAND. The total number of simulated events is 5 x 10® 
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FIG. 8: Angular distribution of final muons at CHOOZ. Tlie total number of simulated events is 5 x 10®. 
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FIG. 9: Comparisons of the cosS and (j) distributions of final muons at KamLAND. The total number of simulated events of 
the present MUSIC simulation is 5 x 10®. The muon fitter results ^9\ represent actual experimental data. The modified MUSUN 
simulation 40j is a standard rock calculation while the present simulation is an exact calculation based on the Ikenoyama rock 
composition. Both the muon fitter and modified MUSUN results are rescaled from arbitrary units to fit the present results. 
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FIG. 



10: Comparisons of the 9 and 

l6 



distributions of final muons at CHOOZ. The total number of simulated events of the 
present MUSIC simulation is 5 x 10". "Simulated Total Flux" refers to the simulated muon flux integrated over the entire 
hemispherical solid angle. The total flux is binned to generate the angular distributions. "Simulated Experiment" is the 
simulated muon flux integrated over a limited range of solid angle described by the geometry of the RFC plates used in the 
cosmic ray experimental setup on the CHOOZ site. Both the experimental and simulated experimental results are rescaled 
such that the small angle parts of all three 6 distributions agree. 
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FIG. 11: Average muon energy E^ versus 8 and (j) at KamLAND. The total number of simulated events is 5 x 10 
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FIG. 12: Average muon energy E^^ versus 6 and (j> at CHOOZ. The total number of simulated events is 5 x 10® 
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FIG. 13: A sketch of a vertical cylindrical detector. The inner volume is indicated by dotted lines which is taken to be 
cylindrical for Super-K and CHOOZ but spherical for KamLAND. 
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FIG. 14; An illustration ol the binning strategy ol slant depth X in the 6 — (j) space. The bars represents regions of the solid 
angle corresponding to the edges of a 3D topographical map and is blocked from the random generation of 6 and 0. The black 
dots represent the original lattice sites from a latticized hill profile. The dotted lines partition the remaining solid angle into 
regions of nearest neighbors. 
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FIG. 15: Integrated muon flux versus muon energy on ground level. The angular integration is talten over the the entire 
hemisphere. Experimental data support the feature that the energy spectrum for _E,i < 1 GeV is almost flat p\. 
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FIG. 16: Integrated muon flux versus the zenith angle on ground level. The energy integration is taken over the range 
0.106-10'' GeV. The experimental data are talten from References |4Ml4^. 
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TABLE I: Integrated muon intensity J^^ and energy i?^ versus vertical depth h from a flat surface in standard rock. Results 
labeled "Sheffield" are taken from Reference p^ that uses the original Gaisser parameterization and MUSIC. Results labeled 
"KSU" present this work using the modified Gaisser parameterization in Eqs. 0-0 and MUSIC. The initial muon energy for 
vertical depth (300 <h< 2000) mwe is (0.106 < S^o < lO") GeV and for (2000 <h< 10000) mwe is (0.106 < E'^o < lO'^) GeV. 
The number of simulated events is lO''. 

Sheffield KSU Sheffield KSU 

h (mwe) J^ (cm-^s-^) J^ (cm'^s'^) g^ (GeV) E^, (GeV) 

500 1.70 X 10'^ 1.71 X 10^5 97 97 

1000 2.20 X lO"*^ 2.21 x 10-"^ 157 158 

2000 1.81 X 10"^ 1.81 X 10^'' 236 236 

3000 2.94 X 10"® 2.95 x 10~** 285 284 

4000 6.33 X 10"^ 6.34 x 10"^ 316 313 

5000 1.58 X 10"^ 1.57 x 10"^ 337 339 

6000 4.30 X 10"^° 4.21 x 10"^° 351 345 

7000 1.24 X 10"^° 1.26 x 10"^° 361 365 

8000 3.73 x 10"" 3.61 x 10"" 369 356 

9000 1.15 X 10~" 1.14 x 10"" 375 373 

10000 3.65 X 10"^^ 3.61 x 10"^^ 380 363 
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TABLE II: Chemical composition of the Ikenoyama and Ardennes rock in elemental percentage. The Ardennes rock composition 
is the average of several samples. The CHOOZ rock data are approximate values only. Details are documented in an internal 
note 35]. 

Ikenoyama Ardennes 

Chemical formula % % 

Si02 60.70 58 

Ti02 0.31 

Ala O3 17.39 19 

Fez O3 1.10 

Fe O 1.22 17 

MnO 0.15 

Mg O 0.93 4 

CaO 6.00 

Na2 O 6.42 

K2 O 3.47 2 

P2 O5 0.18 

H2 O 0.97 

S 0.01 

CO2 0.96 
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TABLE III: Average muon energy E^, muon flux J^, the muon rate inside the tank _R^ and the muon rate inside the inner 
detector volume i?^ for Super-K, KamLAND and CHOOZ. The inner detector of Super-K is a cylinder and that of KamLAND 
is a balloon. The muon rate inside the CHOOZ inner detector is not simulated. 

Site E^ (GcV) J^. (cm~^s~') Rj, (Hz) i^^ (Hz) 

Super-K 271 ± 2 (1.48 ± 0.04) x 10"'^ 2.438 ± 0.004 1.828 ± 0.003 

KamLAND 268 ± 2 (1.70 ± 0.05) x lO"'^ 0.676 ± 0.001 0.246 ± 0.001 

CHOOZ 60.6 ±0.4 (6.12 ± 0.07) x lO"'' 30.5 ± 0.2 
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TABLE IV: Comparisons of muon flux and average energy on the ground level. The experimental values of muon flux is taken 
from Reference |27j . The momentum cut-off of the muon flux measurement is 0.35 GeV. The low total energy cut-off of the 
present calculations is 0.106 GeV. The quoted experimental value [g of vertical muon energy is 4 GeV. The simulated value is 
4.19 GeV. 

Method J^ (cm~^s~^) E^, (GeV) 

Monte Carlo 1.91 x 10"^ 6.92 

Gaussian Quadrature 1.90 x lO"'^ 6.95 

Experiments (1.90 ± 0.12) x 10"^ 



